Thermodynamically self-consistent liquid state theories for systems with bounded 
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The mean spherical approximation (MSA) can be solved semi-analytically for the Gaussian core 
model (GCM) and yields - rather surprisingly - exactly the same expressions for the energy and the 
virial equations. Taking advantage of this semi-analytical framework, we apply the concept of the 
self-consistent Ornstein-Zernike approximation (SCOZA) to the GCM: a state-dependent function 
K is introduced in the MSA closure relation which is determined to enforce thermodynamic con- 
sistency between the compressibility route and either the virial or energy route. Utilizing standard 
thermodynamic relations this leads to two different differential equations for the function K that 
have to be solved numerically. Generalizing our concept we propose an integro-differential-equation 
based formulation of the SCOZA which, although requiring a fully numerical solution, has the ad- 
vantage that it is no longer restricted to the availability of an analytic solution for a particular 
system. Rather it can be used for an arbitrary potential and even in combination with other closure 
relations, such as a modification of the hypernetted chain approximation. 
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I. INTRODUCTION 

It is meanwhile well-known and widely documented 
that conventional integral equation theories - such as 
the Percus-Yevick (PY), the hypernetted chain (HNC), 
or the mean spherical (MSA) approximation - are ther- 
modynamically inconsistent, which means that the vari- 
ous thermodynamic routes to calculate the dimcnsionlcss 
equation of state lead to significantly different results^. 
Over the past decades considerable effort has been de- 
voted to the formulation of thermodynamically self- 
consistent liquid state theories, which, in turn, have led 
to an improved description of the structural and ther- 
modynamic properties of liquids with harshly repulsive 
potentials. In the first generation of these concepts, such 
as the Rogers- Young (RY) 3 , the modified hypernetted 
chain (MHNC)^, or the Zerah-Hansen approach 5 , simple 
functions were introduced in the respective closure rela- 
tions to the Ornstein-Zernike (OZ) equation which use a 
pointwise adjustable but not explicitly state-dependent 
parameter to interpolate between the conventional clo- 
sures. Since self-consistency was enforced for each state 
point independent of the neighboring ones, we shall call 
this approach locally self- consistent. The concepts of the 
second generation of the self-consistent liquid state theo- 
ries were based on more sophisticated ideas: the SCOZA 
scheme^ introduced an explicitly state- dependent func- 
tion in the MSA closure relation to the OZ-equation 
in order to enforce thermodynamic self-consistency be- 
tween different thermodynamic routes; the hierarchical 
reference theory (HRT)£, on the other hand, success- 
fully merged ideas of microscopic liquid state theory and 
renormalization group concepts. In both these advanced 
liquid state approaches thermodynamic consistency was 
enforced in the entire space of system parameters, which 



we shall call global self-consistency. 

In recent years, increasing effort has been devoted 
to investigations of the structural and thermodynamic 
properties of soft matter systems^. The interactions in 
such systems either diverge weakly or even remain finite 
("bounded") at short distances, i.e., when particles over- 
lap. These potentials are commonly referred to as soft 
potentials. Initially, they were investigated by means of 
conventional^ and, more recently, by means of locally self- 
consistent integral-equation theorie aPd? . However, the 
advanced liquid state theories mentioned above have not 
been generalized to systems with soft potentials so far. 
The HRT concept, for instance, relies on the known prop- 
erties of a suitable reference system. While for systems 
with strongly repulsive interactions the hard-core (HC) 
liquid represents an obvious and very successful choice, 
no such reference system can be identified for liquids with 
soft potentials. We therefore have to rule out HRT, at 
least for the time being. 

On the other hand, applications of the SCOZA- 
conceplpi to liquid systems were up to now restricted 
to those cases where the respective interactions can 
be expressed as a combination of HC potentials with 
an adjacent linear combination of Yukawa-tails (HCY- 
systems)**. This restriction can be traced back to the 
fact that the rather elaborate SCOZA-formalism 6 is in- 
tricately linked to the availability of the analytic solution 
of the MSA for such a system^. From this point of view 
the obvious counterpart of HCY systems in soft matter is 
the Gaussian core model (GCM)^ 3 ". For this system the 
pair potential is given by 



$(r)=eexp[-(r/(7) 2 ], (1) 
where e is an energy- and a a length-parameter. 
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Within the framework of the MSA (in the case of 
soft potentials sometimes also termed random phase ap- 
proximation), the structural and thermodynamic prop- 
erties of the GCM can to a large extent be expressed 
semi-analyticallji2*i&. In this contribution we extend the 
SCOZA formalism to the GCM. 



The GCM can be interpreted as a simple model to 
describe soft matter. It has been pointed out that, for 
example, the effective (coarse-grained) interaction be- 
tween two isolated non-intersecting polymer chains or 
dendritic macromolecules can, to a very good approxima- 
tion, be represented by the GCM potentiali^i^. It is for 
two reasons that this potential represents an ideal candi- 
date to apply the concept of global thermodynamic self- 
consistency to systems with soft potentials: (i) as we will 
show, two of the three traditional thermodynamic routes, 
i.e., the energy and the virial route, coincide exactly for 
the GCM within the MSA, a fact that, to the best of our 
knowledge, has not been documented in literature so far; 
there is even evidence^ that this also holds true for other 
systems with soft potentials. Therefore, thermodynamic 
self-consistency has to be enforced between two routes 
only, which considerably facilitates the theory, (ii) Us- 
ing the analytic expressions given by the MSA to enforce 
thermodynamic consistency for the GCM it is possible 
to derive, via standard thermodynamic relations, either 
an ordinary (ODE) or a partial (PDE) differential equa- 
tion for the state-dependent function introduced in the 
closure relation of the SCOZAi£. The ODE results from 
combining the virial and compressibility route and can be 
solved for each isothermal line independently. The PDE 
enforces consistency between the energy and compress- 
ibility route and relates both density- and temperature- 
derivatives. These two differential equations, although of 
different complexity, have to be solved numerically and 
lead within numerical accuracy to consistent results. In 
addition, we propose an equivalent, integro-differential- 
equation (IDE) based formulation of the SCOZA, which 
also has to be solved numerically. This latter approach 
has the advantage that it can be used for an arbitrary 
(soft) potential and in combination with closure relations 
other than the MSA, such as e.g. a HNC-based SCOZA 
ansatz; thus, it is no longer restricted to systems where 
semi-analytic solutions of liquid state theories are known. 



II. MSA 

For the GCM, semi-analytic expressions for the static 
and thermodynamic properties can be derived within the 
MS A 9 ! 16 . Here we add a few details that have not been 
documented yet. 

The MSA closure relation to the Ornstein-Zernike 
(OZ) equation, 



h(r) = c(r) + Q / dr' h{\r - r'|) c(r'), 



(2) 



where h{r) and c(r) are the total and the direct correla- 
tion functions and g is the number density of the system, 
was originally proposed for systems for which the pair 
potential consists of a HC interaction with diameter a 
plus a tail that can take different functional formsi. For 
such potentials, the MSA consists of an ansatz for c(r), 



c(r) = — /3<f>(r) for r > a, 



(3) 



where {3 = (fceT 1 )" 1 , T is the temperature and fee Boltz- 
mann's constant, along with the so-called core condition 
that expresses the impenetrability of the particles 



g(r) = for r < a. 



(4) 



Here, g(r) = h(r) + 1 is the radial distribution function 
(RDF). 

As soft potentials lack a hard core, Eq. cannot be 
applied anymore and the MSA reduces to 



c(r) = — /3$(r) for all r. 



(5) 



For the specific case of the GCM, where $(r) is a sim- 
ple Gaussian, this immediately leads to an analytic ex- 
pression for the static structure factor, S(q), 



S(q) = [1 - gc(q)}- 1 



l + aexp[-(g 2 a 2 /4)]' 



(6) 



where the hat denotes a Fourier transform, q is the wave 
vector, and a = 7r 3//2 p<7 3 /?£. For the RDF we find 



The rest of the paper is organized as follows: in Sec.lTTI 
we re-visit the MSA for the GCM, providing thus the 
basis for the (semi-)analytic formulation of the SCOZA. 
In Sec. IIIII we present the ideas of SCOZA and derive 
the two differential equations and the IDE that impose 
self-consistency and in Sec. IIVI we present details about 
the numerical solution strategies. Sec. |V] is devoted to 
a detailed discussion of the SCOZA results and a com- 
parison with MC simulation data. Finally, in Sec. IVII we 
summarize and draw our conclusions. 



g(r) = 1 - SLJ- [ dqe-** ) ; (7) 
in particular 

ff (0) = l + ^Li 3/2 (-a). (8) 

Here, Li n (a;) is the polylogarithm of order n which is 
discussed in detail in the Appendix. 
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Further, the thermodynamic properties of the GCM 
can be calculated semi-analytically using one of the usual 
three thermodynamic routes^. The results for the dimen- 
sionless equation of state, j3P/ g, where P is the pressure, 
obtained via the compressibility route ('C') 



(9) 



and the virial route ('V') 



where 



H(a) 



Q J 



2a 



a — /3eH(a), 



(10) 



[Li 3 /2(-o0 - Lii/2(-a)] (H) 



have already been reported in2*i&. 

The energy route ('E') has not been considered in the 
literature so far. To obtain (f3P/ g) E we first calculate 
the excess (over ideal gas) internal energy per particle, 

U c */N, 



(3U c: 
~ N~ 



2,t i L , I dr<P(r)g(r)r 2 = |-g [a + Li 3/2 (-a)] 



from which we obtain the excess free energy per particle, 

F cx /N, 



f \J cx ((3',g) a Be r , 



and, finally, the equation of state 



(13) 



III. SCOZA 

The original formulation of the SCOZA for HC 
systems^ is based on the MSA. It enforces the RDF to 
vanish inside the core and for distances larger than the 
core diameter sets the direct correlation function propor- 
tional to the potential; the proportionality factor con- 
tains a state-dependent function that imposes thermo- 
dynamic consistency. Following the same scheme used 
to generalize the MSA to soft potentials (cf. Sec. CJl, wc 
modify the original SCOZA ansatz: 



c(r) = l3K(g, B) $(r) for all r, 



(15) 



where K{g,f3) is an as yet undetermined, state- 
dependent function. As the MSA is recovered for 
K(g,3) = —1, the present formulation of the SCOZA 
takes advantage of the availability of the semi-analytic 
solution of the MSA for the GCM presented in the pre- 
ceding subsection. 

Thus, closed expressions can be derived for the ther- 
modynamic properties within the SCOZA. To sim- 
plify the notation we introduce a function a(g, B) — 
TT 3 / 2 ga 3 f3eK (g, 0) = aK(g,(3), which is explicitly state- 
dependent, but for simplicity suppress the arguments of 
a in the following. 

According to the compressibility route the density 
derivative of the equation of state is given by 



(xL) 1 = (ekBTx$) 1 = l-^(0) = l-a, (16) 

where Xt 1S the isothermal compressibility and the re- 
duced isothermal compressibility xfd ^ s the zero wave- 
vector value of the structure factor S(q). 

Further, following the virial route one finds the follow- 
ing expression for the dimensionless equation of state 



0P 
Q 



1 + 0- 



d (6F C 



dg \ N 



l + -a-0eH(a). (14) 



Thus wc find that virial I jlOp and energy route Ijl4|l 
lead exactly to the same expressions for the dimensionless 
equation of state. This is certainly an unexpected and 
atypical result. In fact, in our numerical investigations of 
similar bounded systems in combination with other clo- 
sure relations, we have observed an analogous, remark- 
able coincidence of the virial and the energy routed. To 
what extent this behavior is a general feature of soft sys- 
tems remains to be investigated. For those systems where 
virial and energy route do coincide, this greatly facilitates 
the formulation of thcrmodynamically self-consistent in- 
tegral equation theories, since consistency has to be en- 
forced only either between the virial and the compress- 
ibility or between the energy and the compressibility 
route. 



PPY 2tt f , 3 d/3$(r) 

— J = l-—Q &rr A ^ g(r) 

o 

1 Be 

= 1+ 2 a ~2~a [ Li5 / 2 ^^ ~ Li 3/2(«)] ( 17 ) 



Finally, according to the energy route the dimensionless 
excess energy per particle is given by 



f3U c: 
~ N~ 



2 a+ 2& L Ll3 / 2 W _a J 



(18) 



The energy and the virial route already coincide within 
the MSA and this also holds for the SCOZA. We are 
therefore left with one single inconsistency, which can be 
removed either via the virial/compressibility or via the 
energy /compressibility route; both possibilities will be 
considered in the following subsections. 
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A. Virial- and compressibility-route 



Eq. ltT7|) with respect to g, 



We start by calculating the compressibility via the 
virial route, Xri which is achieved by differentiating 



J 



(d/3Py _ 1 ( K(g,p) 

V dg J 2ct 3 7t 3 / 2 A>,/3) 2 \ g 



Li 3 /2(a) - Lii/2(o-)] 



dK( e ,p) 
dg 



[2Li 3/2 (5) - Li 5/2 (a) - Li 1/2 (a)] 



(19) 



Equating this result with the compressibility as obtained ODE for the unknown function K(g, (3): 
via the compressibility route (|16fl leads to the following 



dKjg, (3) = K{g, 0) {2ir 3 f3sg 2 a e K(g, g) [K(g, /?) + !]- [Li 3/2 (a) - Li 1/2 (a)] } 
dg g [2Li 3/2 (<5) - Li 5/2 (a) - Li 1/2 (a)] 



(20) 



r 



Note that this ODE can be solved for each isothermal 
line separately. 

Analyzing the ODE, we note that the right hand side 
(RHS) of Eq. (121 III contains two singularities. Obviously 
the denominator vanishes for g — > 0, but expanding nu- 
merator and denominator around g = 0, we find that 



K(g = 0;(3) = 



4V2 



4V2 + f3e' 



(21) 



Further, the denominator also vanishes at a = ao ~ 
—7.7982. This, however, turns out to be a removable 
singularity which can be treated by appropriate means 
(cf. subsection II V A|) . 



B. Energy- and compressibility-route 

To enforce thermodynamic consistency between the en- 
ergy and compressibility route we utilize the variant of 
the SCOZA-formalism proposed iniiLSii which brought 
along a breakthrough of this concept for systems with 
repulsive potentials. This approach is based on replacing 
the differential equation for K (g, (3) by one for the excess 
energy density u — U°*/V. To this end, we consider the 
following thermodynamic relation (see, e.g.jfii) 



Expressing at constant density xj^ d as a function of u, 
the left hand side can be rewritten as 







d 



d(3 Vxf c »y du \ X f cd (u)J 90 
so that finally Eq. 112211 becomes 



^ (23) 



du 
d/3 



d_ 

du 



d 2 u 
dg 2 ' 



(24) 



In contrast to Eq. H20|) . this relation contains deriva- 
tives with respect to both g and f3 and is a PDE of 

the diffusion type. However, the diffusivity, D(g, (3) — 
\ -i -l 

x huj) ' ^ s state-dependent and turns out 
to be negative which renders the numerical solution ex- 
tremely intricate. Xred(u) is now identified with the ex- 
pression obtained by the compressibility route (|lu|) 



du K ^ j 



[x?M] 1 = l-a, 



(25) 



where K(u) is determined by inverting the result of the 
energy route 



d_ 

dp 



Xred 



'dg 2 ' 



(22) 



(26) 
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C. Integro-differential equation approach 

So far, in deducing the SCOZA-ODE JU and PDE 
(|24p. we have taken advantage of the availability of the 
semi- analytic framework provided by the MSA for the 
properties of the GCM. Unfortunately, this represents 
a rather singular exception. In order to eliminate the 
restrictions resulting from this fact one may ask whether 
the SCOZA-concept may be formulated for the general 
case, i.e., when a semi-analytic solution to the MSA is 
not at hand. This is indeed possible as we show in the 
following: let us assume a SCOZA-type closure relation, 



i.e., 

c(r) = fiK§{r) for all r. (27) 

Once K is specified, this leads in combination with the 
OZ equation directly to the radial distribution function 
g(r) = g(r; g, (3; K), which is thus also a function of K. If 
we assume K to be explicitly state-dependent, i.e., K — 
K(g,(3), the compressibility as determined by the virial 
route, i.e., differentiating the standard virial equation of 
state, is 



[gk B Tx) 



. 4tt f 3 d/3$(r) - 2tt 2 . , 



3 d/3$(r) dg(r;g,P;K) 



dr 



do 



2tt 2 dK f 3 d/3$(r) dg(r;g,f3;K) 



3 ^ dg 



dr 



dK 



(28) 



Thermodynamic self-consistency between the virial and 
the compressibility route is now enforced by choosing the 
mixing parameter K such that Xt is equal to Xx> i- e -i by 
finding at fixed temperature T a root of the function 

f(K) = x% - xl- (29) 

Here, derivatives with respect to g and K have to be 
calculated numerically (see below). 

It was exactly this idea that was realized in previous 
applications of parameterized closure relations such as 
RY£, HMSA£, or MHNCl There, however, consistency 
was then achieved only locally, i.e., considering each state 
point in isolation and neglecting thus the density depen- 
dence of K. This corresponds to setting dK/dg = and 
dropping the last term in Eq. I|28|) . In the present ap- 
proach, in contrast, since we consider K to be explicitly 
state-dependent, this term is retained. Thus, the consis- 
tency criterion involves not only isolated state points but 
also via the density derivative nearby state points; there- 
fore we have chosen to call the criterion a global one. The 
quantitative difference between the local and the global 
approaches will be discussed below. 

D. HNC-based SCOZA 

The above approach to thermodynamic consistency 
represents an IDE based re-formulation of the SCOZA; 
it is not only entirely independent of the semi-analytic 
solution provided by the MSA for the GCM and thus be- 
comes completely general in the sense that in this formu- 
lation self-consistency can now be enforced for systems 



with arbitrary (soft) potentials and in combination with 
other, parameterized closure relations. 

To demonstrate the power of this idea we intro- 
duce a HNC-based SCOZA (for clarity we will refer to 
the SCOZA approach introduced above "conventional 
SCOZA" ). This particular choice is motivated by the fact 
that the HNC has been found to work very well for the 
GCM and other soft potentials2ii&. For the closure rela- 
tion of our HNC-based SCOZA we propose 

g(r) = exp(/3^HNc(e,/3)$(r) + h(r) - c(r)), (30) 

where the unknown, state-dependent function 
^hnc(Pj/3) is determined such as to make the RHS 
of the consistency requirement (|29|l . i.e., the equality 
between the compressibility and the virial route, vanish. 
Numerically, this is achieved by solving Eq. I]28p . where 
g(r; g, f3; i^HNc) is now obtained from the solution of the 
OZ-equation along with the closure relation l|30|l . 

IV. NUMERICAL SOLUTION OF THE SCOZA 
AND COMPUTER SIMULATIONS 

A. ODE approach 

The SCOZA-ODE (gDJ has the attractive feature that 
it can be solved for each isothermal line independently 
and represents the fastest route among the three alterna- 
tive formulations presented above to determine K(g,f3). 
We have used an implicit fourth-order Runge-Kutta 
algorithm 2 -, to solve this ODE numerically. This works 
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generally very well except for those state points where 
the expression in the square brackets of the denomina- 
tor of the RHS of Eq. (|2f)|> vanishes. A closer analysis 
shows that this singularity at a = do is removable, since 
both numerator and denominator vanish simultaneously. 
In fact, splitting the density range in two regions, de- 
pending on whether a is smaller or larger than do, and 
integrating the ODE "forward" [starting at g = with 
initial value JHJ] in the former and "backward" (from a 
sufficiently high density so that K = —1) in the latter, we 
were able to smoothly join the partial solutions at d = d Q 
and thus to obtain K (g, (3) over the entire density range. 
We point out that reliable solutions of this ODE can only 
be obtained if an efficient and accurate evaluation of the 
polylogarithm is guaranteed (see Appendix). K(g,(3) as 
a function of g and [3 is displayed in Figure |21 the results 
will be discussed in Sec. Ivl 

Alternatively, we have also solved this ODE with 
MATHEMATICA using a Livermore solver for ordinary dif- 
ferential equations with automatic method switching 
(LSODA) 24 . The poly logarithms encountered in the 
RHS of Eq. ((201) are evaluated in MATHEMATICA with high 
accuracy (for details see^i). Altough the differential- 
equation-solver package is not able to deal properly with 
the removable singularity noted above and breaks down 
for d ~ &o, outside this small range, MATHEMATICA pro- 
vides quasi-exact reference data for the function K (p, /?). 



B. PDE approach 

From the numerical point of view, solving the diffusion- 
type SCOZA-PDE is a delicate task. Boundary con- 
ditions at g — and at a large, but finite p ma x, as well 
as an initial condition at j3 = are required. In par- 
ticular the boundary condition at ^ max has to be chosen 
carefully. In contrast to HC systems, where boundary 
conditions follow naturally from the existence of a max- 
imum density, particles that interact via bounded po- 
tentials can fully overlap and thus it is possible to com- 
press the system to arbitrary high densities. For this 
region, we know that the MSA becomes exact and thus 
self-consistent, i.e., K(g — oo, /?) = —I. In numeri- 
cal calculations, however, we are forced to set K = — 1 
at some finite maximum density g max . Furthermore, we 
have to face the problem of a state-dependent diffusivity 
D(g,/3). Since this quantitiy is even negative, not only 
the solution to the PDE but any numerical error incurred 
in obtaining it may be expected to grow exponentially. 
Among other things, small errors made in the formu- 
lation of the boundary conditions and the inversion of 
the highly non-linear relation l|2t)|) to determine D{g,j3) 
will eventually get dominant. Together, these difficulties 
make it practically impossible to reliably solve this PDE. 
Taking on the other hand K(g,f3) as obtained from the 
ODE and inserting it into Eq. J23J|, we find that 
this relation is fulfilled very accurately, which proves the 
numerical consistency of the two differential equations 



approaches to the SCOZA. 

C. IDE approach 

The IDE based formulation of the SCOZA, i.e., equa- 
tions 129|) along with 128|) . has been solved iteratively 
using both the conventional l|27ll and the HNC-based clo- 
sure (I30|l . We introduce a density-grid (with spacing Ag) 
and assume a starting value K = — 1 . We solve the OZ 
equation with the appropriate closure relation using stan- 
dard integral-equation solver algorithms for a given state 
point (i.e., we fix g and (3) and the neighboring density- 
values (i.e., for g ± Ag). Thus, the derivatives in the 
RHS of Eq. 128|) can be calculated numerically. Due to 
the appearance of the derivative dK/dg, Eq. has to 
be solved iteratively and leads then to K(g,f3) for the 
entire density range considered. As a consequence of the 
iterative and purely numerical character of the solution 
strategy, this approach is more time consuming than the 
solution of the ODE 

D. Monte Carlo simulations 

To test the reliability of our integral equation results 
we have generated reference data for the GCM by means 
of standard Monte Carlo (MC) simulations in the canoni- 
cal ensemble. For each thermodynamic state considered, 
we started from a random configuration of J\f = 1 000 
particles. The system was at first allowed to equilibrate 
for 10 000 passes, where a pass consists of Af trial moves, 
i.e., on average each particle has been subjected to a trial 
move once. After that, we have carried out production 
runs of another 150 — 300 000 passes to calculate the de- 
sired ensemble averages. 

V. RESULTS 

We start the discussion of our results by specifying the 
range in (g, /3)-space where the MSA and the SCOZA 
provide unphysical results, i.e., where g{r) is negative 
(see Figure While this failure of the MSA was briefly 
addressed ini&, we think that a more quantitative analysis 
is in order, since similar problems might be encountered 
in applications of the MSA (and of related concepts) to 
other systems with soft potentials. In fact, also for the 
SCOZA unphysical results can be obtained for certain 
system parameter combinations. For the MSA the limits 
of this range of unphysical behaviour are easily deter- 
mined via Eq. ©, and for the SCOZA they are found 
from the equivalent, generalized expression (i.e., replac- 
ing a by a). Results are shown in Figure ^ indicat- 
ing that at low temperatures the MSA and the SCOZA 
both become unphysical if the density is reduced below 
some threshold density g = g{(3). It is interesting to 
note that similar problems of unphysical solutions and 
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FIG. 1: Region in the density-temperature-space, where the 
MSA and the SCOZA provide unphysical results, i.e., the 
RDF g(r) attains negative values. 




FIG. 2: K(q, /3) as obtained from the solution of the SCOZA- 
PDE 12UH over a representative range of ga 3 and /3e. Note 
that f3e is the inverse reduced temperature, i.e., high values 
of (5e correspond to low temperatures. 

thus restricted applicability have also been reported for 
other self-consistent schemes, such as the RY- or the zero- 
separation concepts, in combination with the GCM 9 . 

The state-dependent function K(g, (3) which guaran- 
tees in the SCOZA-scheme full thermodynamic self- 
consistency, i.e., between all three thermodynamic routes 
is depicted in Figure [21 in a representative part of the pa- 
rameter space. Detailed numerical investigations have 
shown that all three SCOZA formulations presented in 
the previous chapters provide - within numerical accu- 
racy and despite different levels of numerical complexity 
- equivalent results. 

Bearing in mind that the MSA is recovered for 
K(g,(3) = —1, we observe that this function differs sub- 
stantially from this value at low densities (with a pro- 
nounced temperature-dependence) , thus indicating those 
regions where the MSA is thermodynamically inconsis- 
tent. At high densities we confirm earlier results reported 
in2*i&, which have stated that in this regime the MSA 
becomes exact and thus self-consistent. While inpi& this 
conclusion was based on an analysis of the large density- 




FIG. 3: Khnc{q, (3) as obtained from the solution of the 
IDE approach to the HNC-based SCOZA closure 13UH over a 
representative range of ga 3 and (3e. Note that, in an effort 
to enhance the visibility of the data, the viewpoint is now 
different from the one in Figure [5] 

behaviour of the function H as defined in Eq. (|llfl . our 
argumentation follows directly from a visual inspection 
of the function K(g, (3). 

For the HNC-based SCOZA, the corresponding func- 
tion, i?HNC (f?> 0) is shown in Figure^ Taking the devia- 
tion of this function from — 1 as a measure of the thermo- 
dynamic inconsistency of the simple HNC-approach (sim- 
ilar to the case of the MSA), we observe that the HNC 
is to a large degree self-consistent. It is only at small 
densities and low temperatures that -Khnc^j/?) slightly 
deviates from —1. This large degree of thermodynamic 
self-consistency of the HNC for systems with bounded 
potentials was already observed for selected state points 
in2*i£, but was never demonstrated on a quantitative level 
for a wider range of system parameters. 

We conclude our discussion of the thermodynamic self- 
consistency of the conventional (MSA-based) and the 
HNC-based SCOZA concepts by a direct comparison be- 
tween local and global self-consistency, as defined in sub- 
section I^^CJ Let K g (g,(3) denote the explicitly state- 
dependent function K(g,(3), as introduced to enforce 
thermodynamic self-consistency in the IDE formulation 
of the (MSA- or HNC-based) SCOZA (cf. Sec. ImC) : 
thus, the subscript 'g' stands for global self-consistency. 
On the other hand, if the last term in Eq. I|28|) is ne- 
glected thermodynamic self-consistency is only enforced 
for a single, isolated state point and in this case we de- 
note the function by Ki(g,P) (local self-consistency). In 
Figure 0] we show the relative difference between these 
functions for the conventional SCOZA and we observe 
that it amounts to a few percent only for small densities, 
even down to intermediate temperatures. Figure [5] shows 
the same function for the HNC-based SCOZA. Here, the 
differences become noticeable only at small densities and 
low temperatures. Thus, over a large parameter range lo- 
cal consistency is in both cases already a good substitute 
for global consistency. 

We now turn to the structural properties of the GCM 
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FIG. 5: Relative difference between the functions KnNC,g and 
Khnc,i (as denned in the text) over a representative range of 
/3e and ga 3 . Note that, in an effort to enhance the visibility 
of the data, the viewpoint is now different from the one in 
Figure 0] 



by comparing the RDFs for two different thermodynamic 
states. In Figure we have chosen a state-point close 
to the boundary where the SCOZA becomes unphysical 
(cf. Figure We observe that compared to the MC ref- 
erence data, the conventional SCOZA does bring along 
a slight improvement over the MSA. On the other hand, 
the results provided by the HNC and the HNC-based 
SCOZA both reproduce the MC-data perfectly. Figure 
[7| shows the RDF for the GCM at a low temperature 
and low density. Here, we are in the regime where both 
the MSA and the MSA-based SCOZA provide unphys- 
ical results. We see that while the conventional HNC 
results already reproduce the MC data rather well, the 
HNC-based SCOZA leads to a perfect agreement with the 
simulations. We conclude, that although the MSA-based 
SCOZA for the GCM does not bring along the same im- 
provement for the structural properties as documented 
for HCY-systems, the concept of self-consistency by it- 
self proves to be of great value when used with a closure 



FIG. 7: RDF g(r) for the GCM for pe = 10 and ga z = 0.14. 



better adapted to bounded potentials, i.e., a HNC-based 
closure. 

Finally, we examine some thermodynamic properties 
and start our discussion by presenting a rather surprising 
result: if we plot the quantity ua 3 je as a function of the 
density, then the curves evaluated for different isothermal 
lines practically coincide (cf. Figure |8J; even though we 
present only the conventional SCOZA, we note that this 
coincidence is also observed for the MSA, the HNC, and 
the HNC-based SCOZA. This remarkable scaling behav- 
ior might be worth being the subject of future investiga- 
tions. 

We conclude this section with the results for the dimen- 
sionless equation of state, f3P/g, for two different temper- 
atures, i.e., U-qT/e — 10 (see Figure^ and U-qT/e = 0.1 
(see Figures ITU1 and fTTfl . For k B T/e = 10, we find that 
the SCOZA-results coincide with high accuracy with the 
MC-data. For k B T/e = 0.1 we observe (Figure ITU|) that 
the conventional SCOZA provides data that are obvi- 
ously very close to those obtained by simulations, while 
the HNC-based SCOZA data fit them perfectly. A more 
thorough comparison, including this time also other liq- 
uid state theories, such as the MSA, the PY, or the HNC 
approximations^, is displayed on an enlarged scale in 
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FIG. 8: ua 3 /e calculated within the conventional SCOZA as 
a function of ga 3 for different inverse reduced temperatures 
(3e 6 [0, 10]; for discussion see text. 



FIG. 10: (3P/g as a function of ga 3 for the GCM for k B T/e = 
0.1. Note that the conventional SCOZA provides unphysical 
results for the RDF (i.e., g(0) < 0) for ga 3 < 1.27. 
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FIG. 9: f5P/g as a function of ga 3 for the GCM for k B T/e = 
10. The results of the SCOZA and the HNC-based SCOZA 
coincide. Both SCOZA approaches provide physical data for 
the RDF (i.e., g(r) > 0) over the entire density-range. 



FIG. 11: Same as Figure ITHl showing an enlarged view of a 
limited £cr 3 -range. Lines and symbols as labeled. 



VI. CONCLUSION 



Figure ^] We observe that in addition also the virial 
route of the PY and of the HNC (as expectec&iii) nicely 
reproduce the MC data; however, while the SCOZA 
is self-consistent, this is not the case for the conven- 
tional closure relations HNC and PY: their respective 
compressibility-data sometimes differ distinctively from 
their virial and/or energy results. 

Thus we can conclude that the general concept of 
the SCOZA does bring along an improvement over con- 
ventional liquid state theories for a thermodynamically 
consistent description of the properties of the GCM, in 
particular if used in combination with a modification 
of the HNC closure. It is especially remarkable that 
also the structural properties are enhanced, even though 
the SCOZA scheme only enforces self-consistency for the 
thermodynamic properties. 



Motivated by the success of the SCOZA to describe 
the properties of HC systems, we have made first steps 
to extend this concept to systems with soft potentials. 
The fact that the MSA can be solved semi-analytically 
for the GCM makes this system an ideal candidate for a 
first application of the SCOZA. 

Due to the fact that virial and energy route happen 
to yield exactly the same result for the GCM within the 
MSA (and possibly also for other closure relations), we 
are left to fix only one inconsistency, namely between one 
of these two routes and the compressibility route on the 
other hand side. Introducing a state dependent function 
K(g, f3) in the MSA closure, we were able to derive three 
different approaches, namely an ODE, a PDE, and an 
IDE, that enforce thermodynamic self-consistency. While 
the ODE and PDE rely on the analytic solution provided 
by the MSA for this particular system, the IDE formula- 
tion is completely independent of this framework and can 
be applied for arbitrary systems and in combination with 
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any closure relation. It remains to be verified whether the 
IDE approach is also applicable to systems with repulsive 
potentials. 

The three formulations provide results for K(g, (3) and 
K(g, j3) that are equivalent within numerical accuracy. 
In contrast to systems with harshly repulsive potentials 
the improvement of the conventional SCOZA approach 
over the MSA data is less spectacular. While it coincides 
with the MSA results in the limiting case of high densities 
where the MSA is already self-consistent, the (conven- 
tional) SCOZA represents a substantial improvement at 
small densities and low temperatures where the thermo- 
dynamic inconsistency of the MSA is more pronounced. 
Replacing the conventional SCOZA relation by a HNC- 
type closure that contains an analogous state-dependent 
function K-hnq{q, j3), we are able to improve the HNC- 
data for the structural as well as for the thermodynamic 
properties of the system. With this generalised approach 
we have not only demonstrated the flexibility and power 
of the IDE approach but have also proposed what may 
turn out to become a reliable liquid state theory for sys- 
tems with bounded potentials. 

VII. APPENDIX 

The polylogarithm of order n, Li„(z), also known as 
Jonquiere's function 2 ^, is a complex valued function of 
complex argument z, defined by 

oo 

Li„(z) = -— /dt- , (31) 

1 (n) J e l — z 
o 

where n is a positive, real parameter. If z £ K\(l,oo), 
then the polylogarithm is real- valued^. For \z\ < 1 the 
polylogarithm can be evaluated as a power series 



The polylogarithm was introduced irJ^, to calculate the 
thermodynamic properties of the GCM within the MSA 
where, obviously, expression l|32l was used throughout; 
this was done even though there was no guarantee that 
for certain state points the modulus of the respective ar- 
guments \z\ does not exceed 1, violating thus the con- 
dition for the validity of Eq. (|32|) . Since this function 
plays a central role in the formalism of the MSA and 
the SCOZA (see Sec. ITT1 and IIII|) . a reliable evaluation 
of Li„(z) for arbitrary argument z is indispensable for 
a successful solution of the SCOZA-ODE and PDE. We 
therefore provide in the following a more detailed pre- 
sentation of evaluation schemes and indicate how this 
function can be calculated in an accurate and efficient 
way for arbitrary argument z. 

In its evaluation of Li„(z), the MATHEMATICA soft- 
ware relies on Euler-MacLaurin summation, expansions 
in terms of incomplete Gamma functions, and numerical 
quadrature 2 -. Efficient and accurate C- or Fortran-based 
implementations, on the other hand, are more difficult 
to find. First attempts to evaluate Eq. (|31|l directly by 
various numerical integration schemes turned out to be 
either too time-consuming or did not provide results of 
sufficient accuracy. Finally, we found that the following 
functional relation between the polylogarithm and the 
complete Fermi-Dirac function, F n (z), 

oo 

1 f t n 
F n (z) = j d* = -Li„ +1 (-e*) (34) 

o 

along with the accurate and efficient implementation of 
F n (z) via series and asymptotic expansions in combina- 
tion with Chebyshev fits, as implemented in the GNU 
Scientific Library^ 8 , provided the desired results, which 
finally brought the solution of the SCOZA differential 
equations within reach. 



OO 

fc=i 

A relation that turned out to be useful for the present 
application is 

-^-Li„(z) = -Li n -x(z). (33) 
dz z 

A detailed list of additional, helpful relations for this 
function can be found in2£. 



VIII. ACKNOWLEDGMENTS 

The authors are indebted to Elisabeth Scholl- 
Paschinger (Universitat Wien) and Michael Kunzinger 
(Universitat Wien) for useful discussions and to Maria- 
Jose Fernaud for critical reading of the manuscript. This 
work was supported by the Osterreichische Forschungs- 
fonds (FWF) under Project Nos. P15758 and P17823. 
BMM gratefully acknowledges financial support of the 
Erwin Schrodinger Institute for Mathematical Physics 
(Wien), where part of this work was carried out. 



1 J.-P. Hansen and I. R. McDonald, Theory of Simple Liq- 
uids, (Academic, New York, 1986) 2nd edition. 

2 C. Caccamo, Phys. Rep. 274, 1 (1996). 

3 F. J. Rogers and D. A. Young, Phys. Rev. A 30, 999 (1984). 



4 Y. Rosenfeld and N. W. Ashcroft, Phys. Rev. A 20, 1208 
(1979). 

5 J.-P. Hansen and G. Zerah, Phys. Lett. 108, 277 (1985); G. 
Zerah and J.-P. Hansen, J. Chem. Phys. 84, 2336 (1986). 



11 



6 J. S. H0ye and G. Stell, J. Chem. Phys. 67, 439 (1977); 
J.S. H0ye and G. Stell, Mol. Phys. 52, 1071 (1984). 

7 A. Parola and L. Reatto, Adv. Phys. 44, 211 (1995); A. 
Reiner and G. Kahl, Phys. Rev. E 65, 046701 (2002); A. 
Reiner and G. Kahl, J. Chem. Phys. 117, 4925 (2002). 

8 C. N. Likos, Phys. Rep. 348, 267 (2001). 

9 A. Lang, C. N. Likos, M. Watzlawek, and H. Lowen, J. 
Phys. (Condens. Matter) 12, 5087 (2000). 

M.-J. Fernaud, E. Lomba, and L. L. Lee, J. Chem. Phys. 
112, 810 (2000). 

1 E. Scholl-Paschinger, J. Chem. Phys. 120, 11698 (2004). 

2 L. Blum and J. S. H0ye, J. Stat. Phys. 19, 317 (1978); E. 
Arrieta, C. Jedrzejek, and K. N. Marsh, J. Chem. Phys. 
95, 6806 (1991). 

3 F. H. Stillinger, J. Chem. Phys. 65, 3968 (1976). 

4 A. A. Louis, P. G. Bolhuis, J. -P. Hansen, and E. J. Meijer, 
Phys. Rev. Lett. 85, 2522 (2000). 

5 I. O. Gotze, H. M. Harreis, and C. N. Likos, J. Chem. 
Phys. 120, 7761 (2004). 

6 A. A. Louis, P. G. Bolhuis, and J. -P. Hansen, Phys. Rev. 
E 62, 7961 (2000). 

7 B. M. Mladek, M.-J. Fernaud, G. Kahl, and M. Neumann, 
Condens. Matter Phys. 8, 135 (2005). 

8 For the so-called generalized Gaussian core model, 
i.e., $(r) ~ exp[— (r/u)"], n being some arbitrary 
positive number, one can also obtain closed expres- 
sions for selected thermodynamic quantities, which 
are, however, disproportionally more complex than 
for the GCM (i.e., n = 2); for details see B.M. 
Mladek, Integral Equation Theories and Computer 
Simulations for Systems with Bounded Potentials, 
Diploma Thesis (Universitat Wien, 2004, unpublished), 



http : / /tph . tuwien . ac . at/ smt/Mladek/diplomarbeit . pdf 

19 D. Pini, G. Stell, and N. B. Wilding, Mol. Phys. 95, 483 
(1998). 

20 E. Scholl-Paschinger, Phase Behaviour of Sim- 
ple Fluids and Their Mixtures, PhD Thesis 
(Technische Universitat Wien, 2002, unpublished), 
http : //tph. tuwien. ac . at/ smt /paper s/PhD/ESP .pdf. 

21 C. Caccamo, G. Pelicane, D. Costa, D. Pini, and G. Stell, 
Phys. Rev. E 60, 5533 (1999). 

22 Note, that the 'correct' diffusion equation for non-constant 
diffusivity D(g, (3) would read 

du 9 ( t>i a\ du \ 

m = d- e \ D{Q > p) ~d~ e )- 

23 W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. 
P. Flannery, Numerical Recipes in FORTRAN: The Art of 
Scientific Computing (Cambridge University Press, Cam- 
brdige, 1992) 2nd edition. 

24 S. Wolfram, The Mathematica Book (Wolfram Media Inc., 
Champaign, 2003), 5th edition. 

25 Wolfram Research: 

http : / /mathworld . wolfram . com/Polylogarithm . html 

26 D. Wood, The Computation of Polylogarithms,(Technicsl 
Report 15-92, University of Kent, Computing Laboratory, 
Canterbury, 1992). 

27 Wolfram Research: |http : / /functions .wolf ram. com7] 
ZetaFunctionsandPolylogarithms/PolyLog/. 

28 M. Galassi et al, GNU Scientific Library Reference 
Ma-nual (2nd edition), ISBN 09 54161734; see also 
http : //www . gnu. org/ sof tware/gsl 



